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Cosmology and stability in scalar tensor bigravity 
with non-minimal kinetic coupling gravity 

F. Darabf] and M. Mousavf] 

Department of Physics, Azarbaijan Shahid Madani University, Tabriz, 53714-161 Iran 

We generalize the scalar tensor bigravity models to the non-minimal kinetic coupling scalar ten¬ 
sor bigravity models with two scalar fields whose kinetic terms are non-minimally coupled to two 
Einstein tensors constructed by two metrics. We show that a broad class of expanding universes 
can be explained by some solutions of this model. Then, we study the stability issue of the solutions 
by means of imposing homogeneous perturbation on the equations of motion and extract the stable 
solutions. 

PACS numbers: 95.36.+X, 12.10.-g. ll.10.Ef 


I. INTRODUCTION 

The infrared modifications of Einstein general relativity is receiving attention by two aspects: i) the theoretical 
aspect from an effective field theory point of view that is one of the very natural choices to seek when declaring the 
diffeomorphism invariance in GR, and ii) the observational aspect by means of explaining the accelerating expansion 
of the Universe. In this regard, it seems that general relativity has to include a new mass term to validate such 
attempts. Recently, a ghost free 0,1 non linear massive gravity model (the dRGT theory) was constructed M 
involving two metric tensors, one dynamical g llv and one non-dynamical f^ appearing in a unique set of terms. These 
two metrics are used to add a new mass term to the original action (GR), namely the interaction term. For this model 
with a non-dynamical reference metric, it has been shown that the flat Friedmann-Robertson-Walker (FRW) Universe 
dose not exist [H, and that the open FRW solutions arc allowed 0 involving the problems of strong coupling 0 and 
ghostlike instabilities 0. Attempts to screen such problems led people to a natural way of extending the massive 
gravity theory and going beyond it toward a similar new model in which two dynamical symmetric tensors g^ v and 
fitv appear as foreground and background metrics, respectively in a completely symmetric manner [HI [Tb| which is 
called the bigravity theory. Obviously, the massive bigravity theory covers the massive gravity theory and treats with 
two metrics thoroughly in a symmetric way such that convinces one to get rid of the aether-like concept of reference 
metric in massive gravity. Cosmology in massive gravity model has been studied in Refs. EMI, whereas in bigravity 
several branches of regular cosmological solutions have been extracted [20H221. Afterwards, the bigravity model and its 
modification with two independent scalar fields have been investigated in I2.'il-l2~)l| presenting the models which certify 
the stable solution describing the spatially flat FRW solution. In these models, the scalar fields are minimally coupled 
to the metric. As we know, there are varieties of cosmological models which contain scalar fields non-minimally coupled 
to gravity [261428} . Moreover, one can go on the extension of scalar-tensor theories and construct coupling between 
the derivatives of the scalar fields and the curvature, namely the non-minimal kinetic coupling scalar tensor theory 
[29( | . In this work, we consider such non-minimal kinetic coupling scalar tensor bigravity model as a generalization of 
the model studied in [HI . 

The organization of this paper is as follows. In section 2, we review the field equations of motions for the minimal 
bigravity action. In section 3, we construct bigravity models with two scalar fields forming the action of non-minimal 
kinetic coupling scalar tensor bigravity model and again obtain the equations of motions besides extracting the Bianchi 
constraints. In section 4, we show that a broad class of the expansion history of the universe can be explained by 
some solutions of the bigravity model, whereas these solutions do not have always stability against the perturbation. 
In section 5, we go through the stability issue of the solutions. The paper ends with a conclusion. 


* Electronic address: f.darabi@azaruniv.edu 
tElectronic address: mousavi@azaruniv.edu 



2 


II. BIGRAVITY THEORY 


The action of Hassan-Rosen theory named bigravity has the following structure 


Sy,i = Mg / d 4 X\/—detgR < ' 9 ' > + M 2 / d 4 x\/—detfR^' > + 2m 2 M 2 ff / d 4 x\/—detg ^ /3 n e n ( V ff -1 /) ■ (1) 

J J J n=0 

Here g and / M „ are two dynamical tensors in the gravity sector, and R ^ and R^ are the scalar curvatures 
corresponding to the metric tensors g and respectively. It should be noted that M e g is defined as 


1 - 1 1 IT 

m^~mI + m]- U 

The tensor \Jg~ x f is the square root of g" f pv which means, (^J <7 -1 /) p (^J g -1 /) v = g pp f pv = X 9 v . For this 

defined tensor, e„(X)’s are given by 

e 0 (X)=l, ei (X) = [X], e 2 (X) = i ([X] 2 - [X 2 ]) , 

e 3 (X)=i([X] 3 -3[X][X 2 ]+2[X 3 ]), 

e 4 (X) =1 ([X] 4 - 6[X] 2 [X 2 ] + 3[X 2 ] 2 + 8[X][X 3 ] - 6[X 4 ]) , 

efc(X) =0 for k > 4, (3) 

where [A'] means the trace of the tensor X 9 u . For simplicity, we take the minimal and non-trivial case as [3.1] 


Su = Mg J d^x^f—detgR^ + M'j J d 4 x\/—det fR^^ + 2 m 2 M^ s J d 4 x \/—detg ^3 — tn/g 1 / + det\/ff 4 /) , (4) 

which is apparently a direct result of using equation (O in terms of e n as follows 

3 ^ W g~ x f + detV g~ l f = 3e 0 ((Vff -1 /) S) ~ d ((a/^ 31 /) + e4 ((vV 1 /) ^) • (5) 

Considering non-minimal models leads to quite complicated calculations, whereas in the minimal model the interaction 
term of two metrics g pv and f pi/ is just obtained by the trace of (^\Zg~ x f^j w It should be mentioned that this 
simplification does not absolutely change our following results. Before starting the next section, it would be worthwhile 
to explain a little about the details of extracting the field equations of bigravity model. The variation of action © is 
given by 


S g S hi = Mg J d 4 xSg (V-det gR i9) ^ - 2m 2 M 2 ff J d 4 xS g (y/-detg (3 - tr \Jg~ x f + dety/ff -1 /)) , (6) 

and 

<5/Sbi = M 2 J d 4 x5f (^V~detfR^'J — 2m 2 M 2 ff J d 4 xSf det g ^3 — tn/ g~ x f + det• (7) 

Clearly, the first terms of two above equations produce the known Einstein tensors for metrics g pl/ and f pvi respectively 
as 

S g (y/-det = y/-detg (r{?J - ^g pv R (a) ^j , (8) 

and 

S f (V^tetfR^) = v^dit/ (r$ - . 


(9) 
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By considering that <5tr (Vff V) = §tr (s\/s 1 f$9 : ) 
yj det(— g) yj det g~ x f = y / det(—/) we can conclude that 


h tT (j(y/9 1 f) $9 ^ 


beside the property 


0 =Mg (rW - 

+ w 2 M e 2 ff (3 - try's -1 /) + (Vs -1 /) P * + (\/s -1 /) P M |, (10) 

and 

0 =M 2 (r$ - 

+ m 2 M 2 s yJdet (f~ 1 g) |-^/ w (\/3 -1 /) P * - (\/3 -1 /) P m + /^det (Vs -1 /) | ■ (H) 

Here, we are ended up with two independent equations of motion for metrics g ^ and in general model of bigravity. 


III. BIGRAVITY WITH TWO SCALAR FIELD KINETIC TERMS N ON-MINIM ALLY COUPLED TO 

CURVATURES 

Let us review the original scalar tensor bigravity action 


where we have 


=M 2 J d^xyJ-deXgR^ + M] J dfxyj -det. fR U) + 84 , + S^+ 

2 m 2 M 2 s J d 4 xy/-detg (3 - tr \fg~ T j + det y^ -1 /) , 


S(/, = —M 2 J d 4 xy/~detg {egcj)d v <j> + 2 V (<£)} + J d 4 a;£ ma tter(3^, $»), 


( 12 ) 


(13) 


and 


Ss = —M 2 j d 4 xV-det/ {XU.d^d^ + 2 £ 7 ( 0 } , (14) 

where V((f>) and £/(£) are scalar field potentials of metrics g ^ and respectively. Now, we are interested in writing 
the modified form of the above action by means of including the non-minimal kinetic derivative couplings to the 
curvatures in the action as follows 


Stot = M a j d 4 xy/^detg ( R <®> - [eg^ + kG$] - 2V(</>}) + 

Mj J d A x^fctf (i? ( « - [A+ 0G$\ d^d v i - 2 U(0) + 

2 m 2 M 2 s J d 4 xy/-detg (3 - tr \fg~ r ] + detyW) , 


(15) 


where we have ignored the matter contribution, and n and A are the g and / MJ/ coupling parameters with dimension 
of (length) 2 , respectively. Varying the action (fl5l) with respect to and (j) yields the field equations 


0 = - G$ + M 2 T$ + kM 2 &$+ 


2M eS js^ ( 3 - trVff -1 /) + \fw (Vs -1 /) P * + \fv P (t/s -1 /) 


-1 p 


p ( j 


(16) 
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and 


0 = 


'eg^ + nGW 1 


respectively, where V<f, = dV{(j)) / d(f>, V^ 1 = VU and 


7-'(0) — tr 

^ /iy c 


-^(V./*) 2 


- g^V{<p), 


(17) 


(18) 


&jfj = - + 2V a <j>V ( ^R^ a „) + V“0V^i?2 


’[lOLl/fi 


v^v^v.v^- - -G$(V0) 2 + 


9 [tv 


-^v Q v^v a ^v^ + i (n^>) 2 - VaWpR^ 


(19) 


where □ = D ff is the d’Alcmbertian with respect to the metric g. Using (fl5l) . and because / is dynamical as well as 
5 , there is a symmetry between them so that we are allowed to write 


0 = - G$ + M 2 f T$ + /3M o 2 0(/j+ 


^esV det (f~ 1 g) l~Up (\/ff _1 /) P </ - \fu P (vV 1 /^ M + JW det (vV 1 /) 


( 20 ) 


0 = 


Ar+iSG^ 


- t/j, 


( 21 ) 


where L(e = dU(^)/d^. Under the mentioned symmetry condition <-> Vy and □<, f> □/ we just have the changes 


T (?) = A 

/XI/ 


v M ev^--/ M ,(vo 


-UMO, 


( 22 ) 


e$ = - + 2 v^v ( ^i? (/)Q 

V„V“£V„V„£ - V.V^ae - ic^(V0 2 + / M „ i (DO 2 - 


(23) 


According to Bianchi identity we have 


0 = V£ 





(24) 


0 = V? • (25) 

By imposing the covariant derivative Vff and on equations © and and using the constraints 

^MgTffi + nMgQ^^j = 0, and +/3M 2 0^{^ = 0, resulting from equations (fT71) . (fl8l) . (fl9l) . (f2T)l . 

(l22l) and (l23l) . besides refereing to (l24l) and (j25)l . we can obtain two following constraints 

0 = -9^K (irVFv) + l^ g [Up (Vff 31 /)" 1 " „ + Up (v 7 ^ 1 /)’ 1 " „} , (26) 


° = V / 


Vdet(f-ig) | -i (vff 31 /) 1V ^ - 5 (VFU) ^ + /^et (v^/) 


(27) 


which can be used to extract an important constraint on the metric coefficients implying that we can find a dynamical 
cosmology in this model. 
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As an important step to test a theory as a real cosmological model, it would be fruitful to inquire whether it is 
possible to have a model expressing the arbitrary evolution of the expanding universe. As a result, now we consider 
the FRW universe for the g metric besides using the conformal time t = r. We emphasize that in bigravity not 
only the metric gg, v but also /^ is dynamical, therefore we take the following form of metrics 


3 

ds 2 = ^ g^dx^dx” 
0 


a(r ) 2 I —dr 2 + 


,i\2 


i=l 


(28) 


3 3 

ds 2 = ^ f flu dx fl dx l/ = —c(r) 2 dr 2 + b(r) 2 ^^(dx 1 ) 2 . (29) 

fi, u—0 i= 1 

Obviously we are not allowed to take c(r) = 1 nor c(r) = 6 (r), because by this option we will have the Minkowski 
metric for ds 2 that leads to Massive gravity non-dynamical cosmology |32f. Thus, for this case, the (r, r) and (i. i) 
components of (flTH) are given by 


2^2 


0 = -m 2 H 2 - 3m 2 M 2 s (a 2 - ab) + + R(<^)a(r ) 2 M 2 - 9kM 2 ^ 


(30) 


0 =M 2 (2H + iJ 2 ) + m 2 M 2 s (3a 2 - 2 ab - ac ) + 
Hj) 2 3i? 2 </) 2 2H<P4> 


E \- - V^Mt ) 2 I M 2 - 




2 a 2 


Moreover, the (r,r) and (*,i) components of (l20l) yield 


(31) 


2c2 


0 = -3 M}K 2 + m 2 M 2 s c 2 { 1 - ^3 ) + ( + £/(£)c(r ) 2 ) - 9/3M| A ^ 


2 a 2 ’ 


(32) 


0 =M) (-2K - 3AT 2 + 2Ki) + m 2 M 2 ff (-^- + c 2 ) + + £/(£)a(r) 




/3M 2 


£ 2 - 


3AT 2 


+ 3ATL 


V 


2 /f« 


(33) 


by definition of K = b/b and L = c/c. Applying two scalar fields for metrics (l28l) and (|29l) helps us to describe three 
metric coefficients a(r), b(r) and c(r) as three degrees of freedom; this is not possible just by one scalar field. It turns 
out that the equations (l26l) and (1571) carry important results for the variables defined in the metrics (1281) and (l29l) as 
follows 

COu 

cH = bK or — = b. (34) 

a 

This is a constraint imposing on the metrics and relating them to each other. For the case a ^ 0 the above constraint 
gives c = ab/a but for a = 0, we have 6 = 0 meaning that a and b are constant but c is arbitrary. Actually, this 
constraint is allowing us to go on and construct an expanding cosmology. Redefining scalar fields in accordance with 
conformal times g = ( = T as (p = tp(r]) and £ = £(£) with w(r/) = ip'{g) 2 , V(rj) = V {<p{r])), er(£) = £'(£) 2 (prime is 
the derivative with respect to its conformal time) and 17(0 = U (£(£)) facilitate us to rewrite equations (l3Ul) - (l33l) as 
follows 













6 


0 = -3 M 2 H 2 - 3 m 2 M 2 s (a 2 - ab) + 


£UJ ( T ) , 1/2 o.^^Ma) 


+ R(t)«(t) 2 ) M 2 - 9 kM 2 - 


2 a 2 ’ 


(35) 


0 =M 2 (27J + if 2 ) + m 2 M e 2 ff (3a 2 - 2a& - ac) + - A(r)a(r) 2 ^ M 2 - 

[H-3H 2 ^ u(t) + HCj{t) , 


< - 


(36) 


0 = —3M 2 K 2 + m 2 M 2 ff c 2 (l - ^) + + C/(r)c(r) 2 ) M 2 - 9/3M 2 ^^, (37) 

0 =M 2 (-2# - 3 AT 2 + 2ATl) + ro 2 M 2 ff (~^ + c 2 ) + + £/(r)a(r) 2 ^ M|+ 


/3M 2 


3A" 2 \ 

<t(t) ( —A'---h 3KLJ — K&(t) 


Adding and subtracting equations (1551) and (l36l) give us two following equations 


(38) 


0 = 2 M 2 [h — H 2 ^J + m 2 M 2 s (ab — ac) + eco(r)M 2 + 9 w(r) (-6H 2 + H^j + Huj(t) 


(39) 


0 = M 2 ^2H + 477 2 ^ + m 2 M 2 s (6a 2 — 5ab — ac) — 2M 2 a(T) 2 V(r) + [h + 3i7 2 ) w(r) + H6j(t) . (40) 


Again, by subtracting and adding equations (1371) and (1551) we will have 


/ • \ a 3 c a 3 r 2 0Mj r / • \ l 

0 = 2 M) (-AT + KLJ + m 2 M 2 s {~— + —) - Xa(r)M 2 + |ct(t) ( -K - 6K 2 + 3KL J - K&(t) j , (41) 


a 3 a 3 


0 =M 2 (^-6K 2 + 2KL - 2 K) + m 2 M 2 s c 2 ^2 - ^ J + 2 M^c(t) 2 A(t)+ 
/3Mf U(r) (-6AT 2 + 3KL - AT) - Ka{r) 


(42) 


Therefore, for arbitrary metric coefficients a(r), b(r) and c(r), by choosing w(r), a (r), V'(t) and {/(t) satisfying 
equations (1551) - ATI) , we can reconstruct cosmological models with given evolutions of a(r), 6(r) and c(r). 


A. Conformal description of the power expanding universe 


In the last section, in the equations (1251) and (1291) . we have applied the conformal metric with a conformal time r 
to express the evolution of the scalar field 4>(t). In the following, we are going to show how the conformal time can 
describe the known cosmologies. As a result, again we remark the conformally flat FLRW metric as 


ds 2 = o(t) 2 —dr 


2 , V^/j_i\2 


i= 1 


(■ dx l 


(43) 
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Qn 

According to the power expanding universe the scale factor is treating as a(r) = with n / 1. Considering n = 1 
can explain the de Sitter universe in the arbitrary model representing dark energy or inflation. Looking back to the 
power expanding scale factor leads to the following redefinition of the time coordinate as 


so we have 



(44) 


n 

n — 1 

By considering the above discussion, we will have the following form for the FLRW metric 


(45) 


— In. g 

J2(dx 1 ) 2 . (46) 

2=1 

From this equation one finds that for the range 0 < n < 1, the metric describes the Phantom universe j33j . for n > 1 
it describes the quintessence universe and for n < 0 it describes the decelerating universe. For the case of Phantom 
universe, we choose + sign of (l44l) or (14511 and impose the time shift t —>• t — to on (14611 . By this shift, we have the Big 
Rip at t = fo; the present time at t < fo, and the infinite past (t —> —oo) is equivalent to the limit of r —>■ oo. For the 
quintessence range (n > 1) we again choose + sign of (THl) or (H5l) . Hence, the limit of r —>• 0 corresponds to t —> +oo 
and the limit of r —>• +oo corresponds to t —> 0, which may describe the Big Bang. Similarly, for the decelerating 
universe we choose + sign. As a result, the limit of r —>• +oo corresponds to t —>• +oo, and that of r — y 0 corresponds 
to t —^ 0 which again describes the Big Bang. Finally, it should be noted that for the de Sitter case {n = 1) we have 


ds 2 = —dt 2 + 


±(n- 1 .) — 
a 0 


i/r->0 t —> +oo 


and 


if t ±oo t —>• —oo. 


(47) 


B. Dark energy solution with a(r) = b(r) = c(r) 

Our universe can be described by the metric g M „. As a result, we have freedom to choose the functions c(r) and b(r) 
which are not directly describing the expansion of the universe since they are / M „ degrees of freedom in the Einstein 
frame. Any way, we choose the functions c(r) and b(r) equal to a(r) to have more convenient calculation, however, 
we should not deduce that they do not have any physical meaning. Thus, we are allowed to take a(r) = b(r) = c(r) 
leading to K = H = L. As a result, the metric interaction terms in equations (l39l) - (H21) vanish and also by supposing 
Mf = M gi e = A and k = (3 we can obtain 


klo(t) + w(t) 


aQ) 2 e 

H 


k ~6H 



2 a(r) 2 


H 

H 


- H = 0, 


(48) 


i(t) 2 V (t) = 3 H 2 + w(t) 


9 kH' 2 e 


2a(r) 2 2 ) ’ 


(49) 


k<t(t) + a (r) | £ + k |-6 H + ^ 



2«(r) 2 = 0, 


(50) 


and 


a(r) 2 f/(r) = 3 H 2 + a(r) 


9 kH 1 e 


2o(t) 2 2 


(51) 









Obviously, we can deduce that ct(t) = w(r) and so V(t) = U(t). By inserting a(r) = and a(r) = e nr into (H51) 
and (1491) and supposing that ao = l, e = |, k = 1 and cr(r) = w(r) we can find 


w(t) + u>(t) 


f- 


-1 


6n — 1 


V 2nT 2n 1 


2 r 2 


(«■-!) = 0 , 


(52) 


V (r) = 3n 2 r 2 " 2 T w(r) 



r 4n-2 



(53) 


w(t) + w(r) 



- 2 ne 2nT = 0, 


(54) 


V(t) = 3n 2 e _2raT + w(r) Q n 2 e~ 2nT - 1^ . (55) 

Clearly, w(r) and V(t) are indirectly coupled to each other, so the w’s which are obtained from equations (15^1) and 
(l54l) are influenced by the evolution of potential V(t). In order to extract the cosmological constant in our model, 
we need that V (r) become constant. Then, by differentiation both sides of (l53l) with respect to r we can extract the 
expression for w(r) satisfying that of equation (15^1) . Thus, we are facilitated to write w(r) for the case a(r) = and 
V (r) = cte as 


r) = 


(2(1 — n)r 2,l 1 ) ( 


|n 2 r 2n+1 


+ 6 n 2 (n — 1) 


(| n 2 r 2n+l _ 


9n 2 (l—2n)r 4n — 3 +- 


2nr 2n 1 


6n—1 


(56) 


Considering the above result besides referring to (Ifnil) makes it clear that w(t) vanishes and also V(t) becomes 
constant, provided that n = 1. This result is exactly the same as that of scalar tensor bigravity model explaining 
the de Sitter model with a(r) = i as a dark energy universe. In the next section, we will seek the stability of all 
mentioned solutions, even the de Sitter one, under the homogeneous perturbation. 

Now, we consider the scale factor a(r) = e nT . Similar to the previous approach for constant potentials we have 


6n 3 e 1 + 2 ne 2nT (l — |^ 2 e 2nT ) 

(6 n — (| n 2 e~ 2nT — l) + 2 ne _1 — 18n 3 e -1-2rlT ’ 


(57) 


which should take just positive values not to conflict with the definition w(^) = ip'(r]) 2 . Studying the three dimensional 
diagram of the above equation for w(r) with the variables t and n lets us to write 


w(r) >0 —»■ (100 < t < +oo for all values of n). (58) 

Thus, to be more precise, we may describe the cosmological constant solution for the case a(r) = e nT as follows 


for (n > 0 , r —» Too) we have 
for {ti < 0 , t — y Too) we have 


(Vo -T 18n 3 and w(r) —> 4 n 
-T 0 and w(r) -r 0^ . 


(59) 


V. COSMOLOGY AND STABILITY OF SOLUTIONS 

Following the last section, we shall discuss about the validity of these wide range of presented solutions to choose 
just the cases which satisfy both the consistency condition for our model and the stability constraint. In the following, 
we want to investigate the stability of the solutions discussed in the previous section. 
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Going through equations (l48l) and (l50l) . shows us a way to find the positive ranges of oj(t) and <t(t) by means of 
plotting their three dimensional diagrams extracted by equations (H51) and (15U1) in which we put two mentioned scale 
factors a(r) = ^ and a(r) = e nT . As we know, w(r) and ct(t) are p 2 and £ 2 , respectively, so they are not allowed 
to take negative values in order to avoid of ghost and inconsistency in the theory. 

By inserting a(r) = ^ and a(r) = e nT into (H51) and supposing that oo = 1, e = 5 , k = 1 and <t(t) = w(r) we can 
find 

w(r) + w(r) + 2 t 2 ^ 1 (n-l) = 0 , (60) 


6j(t) + oj(t) 


„2 nr 


- 6 n ) - 2ne 2nT = 0 , 


(61) 


respectively. Two above equations have the following solutions 


( r ) in(r)) ci _ 2e |( ( ;! 2 X - 2 ( 6rt ~ 1 ) ln G)) ( n _ 1) ; 

_7 F / 7 r 2 ~ 2Tl \ / r 2-2„ \ 2pfcry 

\^2—2n’ 4n—4n 2 y y n—n' 2 


( 

7 

V 


5n + 2 

2^r 


n 4 (384n 5 + 880n 4 - 500n 3 - 760n 2 - 19n + 15) 

t “2 — 2 n 

e 4 (n-i)n r 2 "-i ( 8 n 3 (48n 3 + 164n 2 + 116n + 15) r 6n + 2n(2n + 5)r 2n+4 + 4n 2 ( 8 n 2 + 26n + 15) r 4n+2 + r 6 ) 


8 n 3 (384n 4 + 1264n 3 + 764n 2 + 4n - 15) 


(62) 


;(r) = 


6nr—• 


4 n 2 C 2 + 2 e' 


6nr —- 




e- 2nr \ 
16n 3 / 



64n 5 



(63) 


where Ci and C 2 are the integration constants. The allowed ranges of r and n resulting in the positive valued ui are 
given by the Table. 1 and Table.2 1 . 


Table 1: Intervals for r and n resulting from the condition ct>(r) > 0 for the case a(r) = -^rr- 


T 

n 

O 

1 

to 

o 

(0.38,0.42) 

(0,1) 

(1, Too) 

(1,10 4 ) 

(0,1) 

(I0 4 ,+oo) 

(0.7, 0.9) 


Table 2: Intervals for r and n resulting from the condition clj(t) > 0 for the case a(r) = e nT . 


1 To find the positive valued lo, we first plotted u;(r, n) in the three dimensional diagrams with respect to r (for three ranges of conformal 
time within two classes of r > 0 and r < 0) and n (for three ranges given in subsection A). Then, we selected the allowed ranges of r 
and n resulting in positive uj. 
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T 

n 

(-1.0) 

(-oo,0.25) 

(-12,-10) 

(-0.57, 0) 

(-0.06,0) 

(1, +oo) 

(0,1) 

(-0.3,0) 

(1,10) 

(0.2,0.73) 

(0,1) 

(0.1, 0.9) 

(0,1) 

(l,+oo) 

(100,+oo) 

(—oo, +oo) 


Referring to equations (l35ll - (1551) and rewriting them by including the conformal times rj and gives us 


0 =2M 2 (h - H 2 ^J + m 2 M 2 s (a(r)6(r) - a(r)c(r)) + eu(r])ri 2 M 2 - 


kMI - 
a(r ) 2 


u>(r])i) 2 ^—6 H 2 + H^j + Hu'(rj)ri 3 + 2Hui(r])firj 


(64) 


0 =M 2 (2H + 4A 2 ) + m 2 M 2 s (6a(r) 2 - 5a(r)6(r) - a(r)c(r)) - 2M g 2 a(r) 2 T/(r ? )+ 


KM 2 r/ • n \ „ , 

9 (H + 3A7 2 ) u}(r])r] 2 + Hui (rj)r] 3 + 2H<jj(r\)r\r\ 


a 


(r) 2 


(65) 


0 =2M ? (-* + KL) 4- (-2^2 + d01!) - MC K’«f + 


/3M 2 

c(r) 2 


ct(C)C 2 (-AT - 6 AT 2 + 3Ai) - ^'(OC 3 - 2Ka(C)CC 


( 66 ) 


0 =Mj (-6K 2 + 2KL - 2K) + m 2 M.V(r) 2 (2 - ^ + 2Af|(7(C)c(r) 2 + 


I3M] 

c(r) 2 


tT(C)C 2 (-6A 2 + 3 A'A - If) - AfT'(C)C 3 - 2Aa(C)CC 


(67) 


Additionally, we can write the scalar field equations as follows 


0 = u}{r])fi + U 4^ + 2Hu>(j])f] + V'{rj)a 2 - ^4 


H 2 oj'(rj)r] 2 


H 2 uj(r))ri + 2HHrju}(ri) 


( 68 ) 


and 


0 = cKCX + + (3^ - £M0C + t/'(C)c 2 - 4 


3 A' 


^'(C)C 2 


cr(C)C) + cr(C)C (9A 3 - 9 k 2 l + 6 kk) 


(69) 


Equations (l48l) and (l49l) will give us the metric coefficients provided that we specify these four equations with a 
function /(r) as 


ee 2 -f( r ) 
,k/'(t?) 


6 /'(f?) 


/'(»?) . 


2(??)e 2/(r) / /"(??) 

« V /'0?) 



= 0 , 




(70) 
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V(r])=3 f' 2 {ri)e 2/(t) + uj(t]) 


§«/»e" 4/(T) 


l p -2/(r) 

2 


(71) 


o-'(C) 


2 e^)/no 
P l /'(C) 


/'(C)) + *(C) 


fAe 2 ^) 

. Pf'iC) 


6/'(C) 


/"(C) 


/'(C) 


= 0 , 


(72) 


C/(C)=3/' 2 (0e- 2/(c) +a;(0 


^/ , 2 (C)e 


-4/(C) 


^ P -mo 

2 


(73) 


As a result, we have 

a(r) = 6 (r) = c(r) = , r) = C = r. (74) 

We intend to study the stability of solution d74l) , so we consider the following perturbation to reconstruct the pertur¬ 
bation equations 


C —> C T <5(), t] —> rj -T < 57 y, b —^ 6 (1 + Sfb) > o —^ a, (1 + 5f a ), A —> K + SK , 77 —> H 5H. 

Having more simplicity persuade us to consider the following terms: 


(75) 


M] = M 2 = = M 2 ; 


(76) 


e = A = 1/2 , k = f 

After a troublesome calculation explained in Appendix A, we obtain 


(77) 


d 

dr 


/ Sr] \ 


( -t 

b 4 

_ &3. 

64 

0 

0 

-h. 

64 

_&2. 

64 

SH 


/l7 | /l6^5 

h2 | hsb3 

0 

h$ 

1 hf$ bi 

h$ | h($ 62 


hi 1 hib4 

hi 1 6164 

hi 

hi hibA 

hi hib4 

SC 


0 

c 4 

C6 

c 5 

Cl 

_£2 

= 

C3 

C3 

C3 

C3 

C3 

SK 


95 

9 i 

96 

92 

54 

53 

Sfa 


0 

1 

0 

0 

0 

0 

\SfJ 


l 0 

0 

0 

1 

0 

0 




( 5ti \ 
SH 

SC 

SK 

Sfa 

\Sf b ) 


(78) 


where the (6 x 6 ) perturbation matrix is called M and its elements have been defined in Appendix A. The eigenvalue 
equation with eigenvalues A has the following form 


0 = A 6 + q 5 A 5 + q A A 4 + q 3 A 3 + q 2 A 2 + qi A + q 0 . (79) 

Since we are interested in finding stable solutions, we should consider the negative value of the eigenvalues A. As 
a result, all the eigenmodes gradually disappear and thus it ends up in a damped perturbation. Note that the 
equation (1791) is respecting the general form of the eigenvalue function expansion of an arbitrary matrix (n x n) as 
/m( A) = (—A) ra + trM(—A )" -1 + ... + detM, which can help us to conclude that this equation with negative eigenvalues 
requires 


<75 = —trM > 0, 


(80) 


or (see the Appendix B) 


trM = - 


h 2 

hi 


Kb 3 

bib A 


h 

b A 


— + 52 < 0 . 

C3 


(81) 
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Similar to the previous argument in finding the allowed ranges of r and n resulting in the positive valued ui, to 
investigate the stable solutions we find the allowed ranges of r and n, given by Table 3 and Table 4, resulting in the 
negative values of trM where use has been made of (1521) and (l63ll with ci = C 2 = 0. Note that, we should be careful 
that the following results are obtained in accordance with the consistency ranges of this model mentioned in Table 1 
and Table 2. For instance, in the case a(r) = ^-we have —trM > 0 in the interval (0.4 < r < 0.6 , 0 < n < 0.1); but 
in this range we face with the negative values of w(r) and hence this interval is ruled out. 


Table 3: Intervals for r and n resulting from the condition —trM > 0 for the case cl(t) — -tt • 


T 

n 

(0.01,0.37) 

[1,2.1) 

(0.42,+ 00 ) 

[1,50) 

(b7) 

(0,0.23) 

(3,7) 

(0.26,0.33) 


Table 4: Intervals for r and n resulting from the condition —trM > 0 for the case a(r) = e nT . 


T 

n 

(0,0.02) 

(-0.28, -0.14) 

(2.8,10) 

(0.2,0.42) 

(0.8,1) 

(0.87,0.9) 

(0,0.35) 

(1,1.3) 

(0.68,0.73) 

(1.01,1.32) 


It would be fruitful to survey the stability of the dark energy solutions (|56l) and (l57l) for the presented scale factors 
a(r) = 7 and a(r) = e nT , respectively. Obviously, the solution (1551) is not stable since the ranges of Table 4 do not 
cover it , thus, we do not have any dark energy solution for the case a(r) = e nT . For the de Sitter solution described 
by (1551) , we should plot trM by using (1551) for the limit n —> 1 + 0 to find the negative values of trM in the following 
intervals for r 


for n —>• 1 + 0 0.13 < r < 0.16 and 1.13 < r < 1.4, 

which is covered with the ranges of Table 3 and so the de Sitter solution of the case a(r) = 7 is stable. 


(82) 


VI. CONCLUSIONS 

In this paper, we have constructed non-minimal kinetic coupling bigravity models with two independent scalar 
fields. It has been shown that a wide range of expansion history of the universe can be explained by a solution of the 
bigravity model, particularly inflation or current accelerating expansion of the universe. This description is different 
from the predictions of the models in massive gravity where the background metric is non-dynamical and does not 
lead to the spatially flat homogeneous FRW cosmological solution. In the original scalar-tensor bigravity theory [23| it 
has been shown that in spite of the difficulty of finding stable cosmological solutions , one can find some explicit stable 
solutions by studying the sign of the trace of perturbation matrix extracted from the perturbed field equations. In this 
work, we have followed the mentioned approach by studying the diagrams of the eigenvalue equation condition (1811) , 
which imposes a constraint on the trace of perturbation matrix M to take just negative values. We considered the de 
Sitter universe in the non-minimal kinetic coupling bigravity model and showed that the current model thoroughly 
endorses the de Sitter universe evolution as well as the scalar tensor bigravity with the choice of a(r) = b[r) = c(r) 
for the solution a(r) = i, and then obtained the cosmological constant V(r) = Vo = A with w(t) = 0. For the 
accelerating solution a(r) = e nT we got the dark energy universe for the ranges obtained in the relation (1551) . 

Also, we considered the constraint —trM > 0 in equation (TSlT) and investigated the stability condition for the ranges 
of r and n mentioned in Table 3 and Table 4. It should be noticed that these intervals are satisfying the mentioned 
allowed ranges for r and n in Table 1 and Table 2. In the first place, we focused on the results of Table 3 for the case 
a(r) = 7 ^-. Seemingly, the first two couple ranges of r and n in Table 3 outline roughly the solution interval r > 0 
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with n > 1 reporting the stability of the quintessence universe as one of the solutions offered in section 4. On the other 
hand, two latter couple ranges of r and n implying the interval 0 < n < 1 in the early universe are interpreting the 
phantom universe. Moreover, we considered the de Sitter solution a(r) = ^ and analyzed its stability in (18211 according 
to the final results reported in Table 3 which surely shows that it is stable because the value n = 1 is included in 
the ranges reported in Table 3 . Furthermore, in the case of other suggested scale factor e nT with the corresponding 
values of r and n mentioned in (1591) . showing the dark energy solution, we cannot distinguish any stability. As a 
result, the solutions of non-minimal kinetic coupling bigravity model has some common stability intervals with that 
of [35j]. Throughout the paper, we just considered the homogeneous perturbation, which is independent of t he sp atial 
coordinates. This is because the inhomogeneous perturbation and/or anisotropic background create ghost [341 ] . and 
the superluminal mode [ 35 J in general leads to the violation of causality [36]. 


Appendix A: The calculations of Eqs. (178l) and (1811) 

In the following, we extract equations Eqs. (l78l) and (I8l1) first by using (l34l) 

£ = ir + f-f- < A1 > 

Having plugged the equations (iTf!) and (IA1D into the equations (1M1) - (1U71) we can eliminate L and c as follows 


0 =2M 2 (H - H 2 ) + m 2 M) s a(r)b(r) (l - + eco^r, 2 M)+ 


a(r ) 2 


;(??)? 7 2 ^—6 H 2 + H^j + Hu;'(ri)r] 3 + 2Hcj(r])f]rj 


(A2) 


0 =2 Mg (h + 2 H 2 ^ + m 2 M) s [ 6a(r) 2 — 5a(r)&(r) — ) — 2M)a(T) 2 V (rj)+ 


g 

kM 2 r 


H 


a 


{rf 


i^H + 3 H 2 ^ uj(r])r) 2 + Hu'(rj)rf + 2 


(A3) 


0 =2 Mi - 


KH 


i{t) 3 K a(r) 3 K 2 


f \ —w + K + 2m M eff ) - 2M0C M f + 


2pM)H 2 
b(r) 2 K 2 


a(C)C 2 ( 2 K + 6AT 2 - 3A'| j - Ka'(() C 3 - 2Ka(()(( 


(A4) 


2 ( o r^2 KH\ 2 2 ( K 2 b(r) 2 cl(t) 3 K 2 

0 =2 M 2 -2K 2 -— ] + m 2 M 2 s [2 w w 


pM 2 

c(t) 2 


a(r) 3 AT 

HJ .' en V H 2 ~ b(r)H 2 ~ b{r)H 

(C)C 2 (-6 K 2 + 3 AT A - k) - Ka'iOC 3 ~ 2tfa(C)CC' . 


2M 2 U(()<t) 


By considering the equations (IA2I) and (IA4I) . we find the following equation 


(AS) 


H _ m2M eS 


r \ I -1 K \ ew ( r ?) ? ? 2 

„ air )b(T) 1- 

2 HM 2 W W \ H 


2 H 2 Ha{r) 


U! 


AT 


2M) 


I _ 1 + K\ Aa(C)C 2 


PH 2 


H J b(r)H 2AT 2 b{r) 2 K 3 

Having eliminated H from (IA2I) and (IA3I) and then from (IA4I) and (IASI) , we calculate 


( 17 ) 1) 2 ^—6 H 2 + H^j + Huj'(j])ri 3 + 2 


u(C)C 2 X ( 2k + 6 AT 2 - 3a| ] - AV(C)C 3 - 2Ka(QCC 


(A6) 
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9 k / 7 2 \ 

Mry) 


a (r) 2 V{ri)^ = 3 H 2 + m ^ eff (3a(r) 2 - 3a(r)6(r)) , 


(A7) 



9/5T7 4 \ 
2b{r) 2 K 2 ) 


b(r) 2 U( C)) 


3/7 2 + 


? ” 2M eff 




Combination of (1A61) . (IA7I) and (IA 8 I) gives us the following equation 


(AS) 


0 ={K-H) 


U(Ob(T) 2 K 
2 H 2 


V(v)a(r) 2 

2 77 


K 

4a(r) 2 




+ u)'{rj)fj 3 + 20 ^( 77)7777 


7 2 M , 2 


eff 


-A / «(t) 3 _ , 2 \ _ 3 / a(r ) 2 _ 2 \ / _ K \ ( a{rfH _ a(r)b(r) \ 

2H 2 M 2 \ b(r) [ ’ ) 2HM 2 \ b(r) / V H / \4&(0 - ff2 4/7A / 2 / 

/2 /2 KH 2 6/7 2 3/7/7 \ a'{OC 3 H 2 2H 2 a{()(( 

CTC ^ _ A3 _ + _ A J 7^ 7^ ' 


P 

46(r) 2 

(A9) 


Regarding the independent equations (IA3I) . (1A7I) . (IA 8 I) and (IA9D beside using (l64l) . (1651) . (l 66 l) and m help us to 
extract the following four perturbation equations 


0 =s„ (2+ ASA + 5H {10H+ + 

\ a(r) z J \ ti a \ r ) “ 2 H Ha(r) z 

2m 2 a(r) 2 ( 2 2 2 , , /_ 36k/7 2 36k/7 3 \ 

y -<5/7 + Sf a 12 ?n 2 a(r ) 2 - 167/ 2 + a; ( 77 ) 2- 7— ^ H-r— 7 - ) + 477 ) - 


77 


12m 2 a(r) 2 5 fb + Si 7 


a(r ) 2 a(r ) 2 / 


x -2/7 - 


uj{i])nH 12ui(ri)KH 2 
k /7 


2H 2 - 27/ - 


-’(0 


At// 


i ( r ) 2 / 7 a ( r ) 2 


a(r ) 2 

12 k //// 36 k // 3 77 


;(t7)(57? 


z(t) 2 


z(r) 2 


2H 


Sri — 14/7 7 / + 16// 3 — 2/7 


2 H 2 

IT 


(A10) 


0 =Sf a ^6Z7 2 + 07 ( 77 ) ^ ^^ - 6m 2 a(r) 2 j + Sf b (6 m 2 a{r) 2 ) + SH 


9nHuj(ri) 

a ( T ) 2 


-67/ 


.. , 1 9 k // 2 

+ ^ l 2^^F 


x 07 ( 77 ) + Srj u)(rj) 


( 9 k // 7 / 18 k // 3 77 


V °( r ) 2 


o(t) 2 


y I + 6/77/ - 6/7 3 


(All) 


0 =<J/ 0 (-6m 2 a(r) 2 ) + <5/ fe ^6/Z 2 - yy + + 6m 2 a(r) 2 ) + (5( x 




27/ 


a(r) 2 


(57/ ( - 


a(r) 

07 ( 77 ) 9ku}(t])H\ 


07(77) 9nu}(r])H 2 


a(r) 2 


27/ 


r V , , / 9 k // 7 / 18 k // 3 H 


6/7/7 - 67/ 3 
(A12) 


0 = (<5/7 - (5/7) 


^w(ry) 


15 k 


+ 


k /7 


4a(r ) 2 2a(r) 2 H 2 87/ 2 


3 7 / \ 

2 + H 2 ) 

(<507(77) - MO) 


9 k /7 

4 a ( r ) 2 
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+ (Sv - S() w(r/) 


(Srj-SC) (u(v) 


3 kH 9 kH 2 


a{Ty 


a(r ) 2 


3k!7 


1 

2 ~~ 
kH 


H 


+ 


kH 


kH 2 


8H 2 4a{r) 2 H 4a{r) 2 H 2 


H__H 
2H 2 


(<Stf - (W) 


H H 2 
2H ~ 2H 2 

kuj(t]) 


7 H 


-4 H 2 


8 H ' 2a(r ) 2 4 a(r) 2 H I ' 2H 2 I ' V“ ““7 2a{r) 2 H' 

Apparently, the equation (1A11I) determines Si) by which we can extract SH from (IA10I) as 

Sv = ^Sf a - fSf b - ^SH - T-Sr). 

04 O4 O4 O4 


(A13) 


(A14) 


SH = SH - 


hi /1164 
Equation (IA12I) gives us <5£ as 


!K 1 ~t) + «■ (-r + rr) + (~r + tt) + s ” (~t 

hi) V hi hibij \ hi hibi) V hi 


<5C = 5fa - —Sfb - —SH - —SK - —6(. 

C3 c 3 c 3 c 3 c 3 


h§b 3 

hi 64 


Finally, in order to obtain SK, we substitute the equations (IA14I) - (IA16I) into (IA13I) as follows 


where 


SK = giSH + g 2 SK + g 3 Sf b + gkfa + g^Sg + 


_ h 2 h 3 b 3 b 3 di 

hi hibi bid 3 

d\ 

d? 

b 2 di 


Cidi di 
c 3 d 3 e ?5 ’ 


9 2 = 


h 3 c 3 di 
hi c 3 d 5 


+ 


c 2 d± 


h 3 h 3 b 2 

hi hibi b/id 3 c 3 d 3 
hi h 3 bi bidi Cidi 

hi hi 64 64^5 c 3 d 3 

_ hr heb 3 b 3 di d 3 

hi hi 64 bid 3 d 3 ’ 


d 2 

<V 

(h 

da’ 


9 6 = 


06^4 d 3 

c 3 d 3 d 3 ’ 


(A15) 


(A16) 


(A17) 


(A18) 


and also 


hi =2 


KOJ 
o ? 


h 2 = 10 H 


2 m 2 a 2 V2nHu> 2 H 


H 


2„2 


h 3 = - 


2 m 2 a 
H 


U! 


kloH 


H 2H a 2 H ’ 


hi = 12 m 2 a 2 — 16 H 2 + w ( 2 — 
he = — 12 m 2 a 2 , 


36 kH 2 


4 H, 


he. — I 


kH 12 kH 2 1 


+ 


2 H 2 - 2 H, 


(A19) 


bi =6 H 2 + u 


2 / 18 kU 2 1 i G 2 2 


— — — 6m a 
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2 2 


62 = 6 m a 

9 kw77 

b 3 = --- 6 77, 


bi = oj ( — — 


, H 

b 5 = w — 


1 9 k77 2 


a‘ 

9 k7777 

~ 


18k77 3 \ 


ci = - b 2 , 
C2 = -&i, 
1 
2 


c 3 = W - - 


9k77 2 


C 4 = W 

C 5 = W 

C6 = b 3 


1 

18k77 

277 

a 2 

-1 

9k77\ 

277 + 

a 2 7 


- 677 , 


(A20) 


(A21) 


d± = ui 


— 15k 


k77 


4a 2 2a 2 77 2 877 2 


d 2 = 477 - J 


/ 9k77 




1 


2a 2 277, 


d 3 =io 


f 3k77 —9k77 2 1 77 

l “a 2 a 2 + 2 ~~ Aff 2 


3 77 

2 + 7T 2 ’ 


k 77 
+ 4a 2 77 


k77 2 \ 
4a 2 77 2 J 


77 

277 ~~ ' 

C?4 = W 


77 2 777 „ rr2 9k77 

27T 2 + Y” _ + 1 Y~ 


3kT7 k77 


2a 2 4a 2 77 


F(a), 

77 
‘ 277 


77 

T’ 


(A22) 


where <5er(C) = 7^(a)(5£, 8u}(r)) = F(a)6r] and 7 r (o) is of the order a 3 or 77 3 . Herein, looking back to the equations 
(H51) and (1501) shows that we cannot work with the analytical exact expression for w(t) and er(r) unless we specify the 
scale factor a(r) by which we become able to solve (B81) . Any way, using the equations (|A14I) - (IA17I) beside the two 
following equations 


<577 = Sf a , SK = Sf b , (A23) 

we seek the range of stability of the solutions. It is worthwhile to note that w(t) appears in the final form of the 
relations (IA14I) . (1A151) . (IA161) . (IA17I) and (IA23D in which we replace it by (l62l) and (l(J3l) (it should be remarked that 
we have chosen ci = C 2 = 0 ) for a(r) = rpr and a(r) = e nT , respectively. 


Appendix B: The expressions for tr M 


For the trM, we have 
trM = (24a 10 (w - 477 2 + 277) + 

72uj 3 H 2 (—477T7 4 + 77 + 377 2 (4 + 2777)) + 2a 8 (l48877 4 + w (l3w + 747?) - 477 2 (l69w + 2587?)) + 
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6 uj 2 o 2 H 2 (-6a; - 11052F 4 + (8 - 57 lo) H + 3 H 2 (l36 + 223a; + 924tf) ) + wa 4 x 

(-26928 H 6 + 5 uj 2 H + 96 H 4 (l6 + 187a; + 114#) - 77 2 (a; (96 + 157a;) + 36 (8 + 79a;) 77) ) + 2a 6 77 2 x 

144077 4 + Jp ( w + 36#) + 124 ^ 2 ( 8 + 839w + 84i ^) - ( 24 + 613a;) + ^pr ( 8 + 223w )j ^ x 

(277a; (w + 2a 2 ) (a 2 - 1877 2 ) (12a 4 - 3677 2 a; + a 2 (w - 1277 2 ))) _1 . (Bl) 
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